Self-Organization of Mobile Populations in Cyclic Competition 



Tobias Reichenbach ^, Mauro Mobilia and Erwin Frey ^ 

^Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience, Department of Physics, 
Ludwig-Maximihans-Universitat Miinchen, Theresienstrasse 37, D-80333 Miinchen, Germany. 
^Mathematics Institute and Warwick Centre for Complexity Science, The University of Warwick, Gibbet Hill Road, Coventry 

CV4 7AL, United Kingdom 



Abstract 



00 

O 

o 

(N 

s 

00 
(N 



The formation of out-of-equilibrium patterns is a characteristic feature of spatially-extended, biodiverse, ecological sys- 
tems. Intriguing examples are provided by cyclic competition of species, as metaphorically described by the 'rock- 
paper-scissors' game. Both experimentally and theoretically, such non-transitive interactions have been found to induce 
self-organization of static individuals into noisy, irregular clusters. However, a profound understanding and characteriza- 
tion of such patterns is still lacking. Here, we theoretically investigate the influence of individuals' mobility on the spatial 
structures emerging in rock-paper-scissors games. We devise a quantitative approach to analyze the spatial patterns 
self-forming in the course of the stochastic time evolution. For a paradigmatic model originally introduced by May and 
Leonard, within an interacting particle approach, we demonstrate that the system's behavior - in the proper continuum 
limit - is aptly captured by a set of stochastic partial differential equations. The system's stochastic dynamics is shown 
to lead to the emergence of entangled rotating spiral waves. While the spirals' wavelength and spreading velocity is 
demonstrated to be accurately predicted by a (deterministic) complex Ginzburg-Landau equation, their entanglement 
results from the inherent stochastic nature of the system. These findings and our methods have important applications 
for understanding the formation of noisy patterns, e.g., in ecological and evolutionary contexts, and are also of relevance 
for the kinetics of (bio)-chemical reactions. 
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Spatial distribution of individuals, as well as their mobil- 
ity, are common features of real ecosystems that often come 
paired (May 19741. On all scales of living organisms, from 
bacteria residing in soil or on Petri dishes, to the largest 
animals living in savannas - like elephants - or in forests, 
populations' habitats are spatially extended and individu- 
als interact locally within their neighborhood. Field studies 
as well as experimental and theoretical investigations have 
shown that the locality of the interactions leads to the self- 
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20061). Another important property of most 
individuals is mobility. For example, bacteria swim and 
tumble, and animals migrate. As motile individuals are ca- 
pable of enlarging their district of residence, mobility may 
be viewed as a mixing, or stirring mechanism which "coun- 
teracts" the locality of spatial interactions. 

The combined influence of these effects, i.e. the com- 
petition between mobility and spatial separation, on the 
spatio-temporal development of populations is one of the 
most interesting and complex problems in theoretical ecol- 
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bach et al. 2007a I. If mobility is low, locally interacting 
populations can exhibit involved spatio-temporal patterns, 
like traveling waves (Igoshin et al. 20041, and for exam- 



ple lead to the self-organization of individuals into spirals 



in myxobacteria aggregation (Igoshin et al. 20041 and in- 



sect host-parasitoid populations (Hassell et al. 1991), or 
more fractal-like structures in competing strains of E.coli 
(Kerr et al. 2002|. On the other hand, high mobility re- 



sults in well-mixed systems where the spatial distribution 



of the populations is irrelevant (Maynard Smith 1982 Hof- 



bauer and Sigmund 1998). In this situation, spatial pat- 
terns do no longer form: The system adopts a spatially 
uniform state, which therefore drastically differs from the 
low-mobility scenario. 

An intriguing motif of the complex competitions in a pop- 
ulation, promoting species diversity, is constituted by three 
subpopulations exhibiting cyclic dominance. This basic mo- 
tif is metaphorically described by the rock-paper-scissors 
game, where rock crushes scissors, scissors cut paper, and 
paper wraps rock. Such non-hierarchical, cyclic competi- 
tions, where each species outperforms another, but is also 
itself outperformed by a remaining one, have been identified 



in different ecosystems like coral reef invertebrates (Jackson 



and Buss , 1975 1, rodents in the high-Arctic tundra in Green- 



land (Gilg et al. 2001 1, lizards in the inner Coast Range of 



California ( Sinervo and Lively 1996 ) and microbial popula- 



tions of colicinogenic E. coli (Kerr et al. 2002 Kirkup and 



Riley 2004). In the latter situation, it has been shown that 
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Figure 1: The stochastic spatial system at difFerent scales. Here, each of the colors yellow, red, and blue (level of gray) represents 
one species, and black dots identify empty spots. Left: Individuals are arranged on a spatial lattice and randomly interact with their 
nearest neighbors. Middle: At the scale of about 1,000 individuals, stochastic effects dominate the system's appearance, although 
domains dominated by different subpopulations can already be detected. Right: About 50,000 mobile interacting individuals 
self-organize into surprisingly regular spiral waves. 



spatial arrangement of quasi-immobile bacteria (because of 
'hard' nutrient or substrate) on a Petri-dish leads to the 
stable coexistence of all three competing bacterial strains, 
with the formation of irregular patterns. In stark contrast, 
when the system is well-mixed, there is spatial homogene- 
ity resulting in the take over of one subpopulation and the 
extinction of the others after a short transient. 

It is worth noting that the emergence of noisy patterns, 
as those studied here, is a feature shared across disciplines 
by many complex systems characterized by their out-of- 
equilibrium nature and nonlinear interactions. Examples 
range from the celebrated Belousov-Zhabotinsky reaction 



(Zaikin and Zhabotinsky 19701 (spiralling patterns) and 



many other chemical reactions (R. Kapral and K. Showal- 
ter 



1995[), to epidemic outbreaks (traveling waves) (Gren- 
2001 Cummings et alT] 2004[), excitable me- 



fell et al 



dia (Muratov et al. 2007 R. Kapral and K. Showalter 



19951, and calcium signalling within single cells (Lechleiter 



et al. 1991 Falcke 2004 Bootmann et al. 2006). More 



over, cyclic dynamics as described by the rock-paper-scissors 
game occur not only in population dynamics but have, e.g., 
been observed in social dilemmas relevant in behavioral sci- 
ences ( Sigmund et al. 2001 Hauert et al. 20021. Therefore, 



we would like to emphasize that the methods presented in 
this work are not limited to theoretical ecology and biology, 
but have a broad range of multidisciplinary applications and 
notably include the above fields. 

Pioneering work on the role of mobility in ecosystems was 



performed by Levin (19741, where the dynamics of a pop- 



ulation residing in two coupled patches was investigated: 
Within a deterministic description. Levin identified a criti- 
cal value for the individuals' mobility between the patches. 
Below the critical threshold, all subpopulations coexisted, 
while only one remained above that value. Later, more re- 
alistic models of many patches, partly spatially arranged, 
were also studied, see Hassell et al. (1991 1994 1; Blasius 



et al. (19991; D. Alonso (20021 as well as references therein 



These works shed light on the formation of patterns, in par- 
ticular traveling waves and spirals. However, patch models 
have been criticized for treating the space in an "implicit" 



manner (i.e. in the form of coupled habitats without inter- 
nal structure) (Durrett and Levin 1998). In addition, the 



above investigations were often restricted to deterministic 
dynamics and thus did not address the spatio-temporal in- 



fiuence of noise. To overcome these limitations, Durrett and 



Levin (1997) proposed to consider interacting particle sys- 



tems, i.e. stochastic spatial models with populations of dis- 
crete individuals distributed on lattices. In this realm, stud- 
ies have mainly focused on numerical simulations and on 
(often heuristic) deterministic reaction-diffusion equations, 
or coupled ma ps (,Durrett and Levin 1994 1997 1998 , King] 



and Hastings 2003| Czaran et al. 2002 Liebermann et al. 
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Here, we demonstrate how a - spatially explicit - stochas- 
tic model of cyclically interacting subpopulations exhibits 
self-formation of spatial structures which, in the presence of 
individuals' mobility, turn into surprisingly regular, geomet- 
ric spiral waves. The latter become visible on the scale of a 
large number of interacting individuals, see Fig. [l] (right). In 
contrast, stochastic effects solely dominate on the scale of 
a few individuals, see Fig. [l] (left), which interact locally 
with their nearest neighbors. Spatial separation of sub- 
populations starts to form on an intermediate scale. Fig. [T] 
(middle), where mobility leads to fuzzy domain boundaries, 
with major contributions of noise. On a larger scale. Fig. [l] 
(right), these fuzzy patterns adopt regular geometric shapes. 
As shown below, the latter are jointly determined by the de- 
terministic dynamics and intrinsic stochastic effects. In the 
following, we elucidate this subtle interplay by mapping - in 
the continuum limit - the stochastic spatial dynamics onto a 
set of stochastic partial differential equations (SPDE) and, 
using tools of dynamical systems (such as normal forms and 
invariant manifolds), by recasting the underlying determin- 
istic kinetics in the form of a complex Ginzburg-Landau 
equation (CGLE). The CGLE allows us to make analyti- 
cal predictions for the spreading velocity and wavelength of 
the emerging spirals waves. Below, we provide a detailed 
description of these methods and convey a thorough discus- 
sion of the spatio-temporal properties of the system with an 
emphasis on the role of spatial degrees of freedom, mobility 
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and internal noise. 



In our first article on this subject (Reichenbach et al. 



2007a) we have described how a mobility threshold sepa- 



rates a biodiverse regime (arising for low mobilities) from a 



high-mobility regime where diversity is rapidly lost . In ( Re- 
[ichenbac het al.[|2007b l we have further analyzed the travel- 
ling spiral waves that arise for low mobilities and computed 
correlation functions as well as the spirals' wavelength and 
spreading velocity. In this article, we provide a compre- 
hensive discussion of the quantitative analysis of the sys- 
tem's properties. This includes the detailed derivation of 
all mathematical equations, an accurate description of the 
numerical simulations (via the implementation of an effi- 
cient algorithm for the lattice simulations taking exchange 
processes into account) as well as the analytical treatment 
of the out-of-equilibrium patterns emerging in the course of 
the time evolution. 



2 Simulations and Results 

We study a stochastic spatially-extended version of a three 
species model originally investigated (on rate equations 
level) by May and Leonard^ (1975| . In ( Reichenbach et al 



2007a ) we have already considered some properties of such 
a model and demonstrated the existence of a critical value 
of the populations' mobility separating a biodiverse regime, 
where all subpopulations coexist, from a uniform regime, 
where only one subpopulation survives. A brief account of 
the analysis of the spatio-temporal properties of the coexis- 



tence phase has recently been given in (Reichenbach et al 



2007b I. Here, we complete and considerably extend those 
previous works by giving a comprehensive picture of the sys- 
tem's properties and details of various mathematical meth- 
ods (interacting particle systems, stochastic processes, dy- 
namical systems, partial differential equations) allowing to 
deal with these kinds of multidisciplinary problems. 

Consider three subpopulations A, B and C which cycli- 
cally dominate each other. An individual of subpopulation 
A outperforms a B individual through "killing" (or "con- 
suming"), symbolized by the ("chemical") reaction AB 
A(Z), where denotes an available empty space. In the same 
way, B outperforms C, and C beats A in turn, closing the 
cycle. We refer to these processes as selection and denote 
the corresponding rate by a. To mimic a finite carrying 
capacity, we allow each subpopulation to reproduce only if 
an empty space is available, as described by the reaction 
ACd AA and analogously for B and C. For all subpop- 
ulations, these reproduction events occur with rate fi, such 
that the three subpopulations equally compete for empty 
space. To summarize, the reactions that define the model 
(selection and reproduction) read 



AB 


-^A0, 


A(Z) 


^AA, 




BC 


^B(Z), 


B(Z) 


-^BB, 




CA 


^C0, 


C(d 


CC. 
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Figure 2: The phase space of the non-spatial system. It is 
spanned by the densities a, b, and c of species A, B, and C. 
On an invariant manifold (in yellow/light gray, see text), the 
flows obtained as solutions of the rate equations Q (an example 
is shown in blue/dark gray) initially in the vicinity of the reac- 
tive fixed point (red/gray, see text) spiral outwards, approaching 
the heteroclinic cycle which connects three trivial fixed points 
(blue/dark gray). In Subsection 3.2 
priate coordinates {yA,yB,yc) 



we introduce the appro- 



which reveal the mathematical 
structure of the manifold and reflect the cyclic symmetry of the 
system. 



these reactions are embodied by rate equations for the tem- 
poral evolution of the mean densities a{t) , b{t) , c{t) of the 
subpopulations A, B and C, respectively, given by Eqs. Q 
in Subsection |3.1[ These equations provide a determinis- 
tic description which is well suited to describe well-mixed 
systems with a large number of individuals, such as Moran 

i2005 2006; ) or urn 



processes (Moran 1958 Traulsen et al 



models (Feller 1968 Reichenbach et al 



20061. For the sys- 



In the absence of spatial degrees of freedom, the kinetics of 



tem under consideration, the rate equations are given and 
carefully investigated in Subsection |3.1[ As main features, 
they possess three absorbing fixed points corresponding to 
survival of only one subpopulation (the solution correspond- 
ing to an empty system is also an absorbing state, but is 
irrelevant for our purpose and will be ignored thereafter), 
as well as a reactive fixed point where all three subpopu- 
lations coexist, see Fig. |2] The coexistence fixed point is 
unstable, and trajectories starting in its vicinity spiral out- 
wards. The spiralling fiows lie on an invariant manifold, 
and approach a heteroclinic cycle which connects the three 
absorbing fixed points. The approaching trajectories spend 
longer and longer periods of time in a neighborhood of each 
(absorbing) fixed points, before departing to the next one. 
This means that the system alternates between states where 
nearly only one of the three species is present, with rapidly 
increasing time period. However, this picture is idealized 
as it crucially (and tacitly) assumes the presence of an in- 
finite population. In fact, fiuctuations, e.g. stemming from 
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a finite number of individuals, lead to reach one of the ab- 
sorbing states where only one subpopulation takes over the 
whole system. Which of the absorbing states is reached 
depends on the initial conditions as well as on random fluc- 
tuations. Due to the symmetry of the reactions ([l]), when 
averaging over all possible initial conditions as well as fluc- 
tuations, all species have an equal chance to survive. 

In the spatially-extended stochastic version of the model, 
we adopt an interacting particle description where individ- 
uals of all subpopulations are arranged on a lattice. In this 
approach, each site of the grid is either occupied by one 
individual or empty, meaning that the system has a flnite 
carrying capacity, and the reactions ([!]) are then only al- 
lowed between nearest neighbors. In addition, we endow 
the individuals with a certain form of mobility. Namely, 
at rate e all individuals can exchange their position with 
a nearest neighbor. With that same rate e, any individual 
can also hop on a neighboring empty site. These exchange 
processes lead to an effective diffusion of the individuals. 
Denote L the linear size of the d-dimensional hypercubic 
lattice (i.e. the number of sites along one edge), such that 
the total number of sites reads N = L"^. Choosing the lin- 
ear dimension of the lattice as the basic length unit, the 
macroscopic diffusion constant D of individuals stemming 
from exchange processes reads 



(2) 



D = edr^N-^''^ : 



the derivation of this relation is detailed in Subsection 13.41 
How do individual's mobility and internal noise, in addi- 
tion to nonlinearity, affect the system's behavior ? 
Insight into this ecologically important issue can be gained 
from a continuum limit where the diffusion constant D is fi- 
nite. Namely, we investigate the limit of infinite system size, 
N oo, where the diffusion D is kept constant (implying 
that the local exchange rate tends to infinity, e —^ oo). In 
Subsections 3.4 and 3.5 we show how. 



in this limit, a de- 
scription of the stochastic lattice system through stochas- 
tic partial differential equations (SPDE) becomes feasible. 
These SPDE describe the time evolution of the (spatially 
dependent) densities a{r,t),b{r,t),c{r,t) of the subpopu- 
lations A, B, and C, respectively. The expression of the 



SPDE, Eqs. (30 1, is given in Subsection 3.5 along with their 



detailed derivation. The latter relies on a system-size ex- 
pansion in the continuum limit which allows to obtain the 
noise terms of the SPDE. Noise stems from the stochasticity 
of the reactions ([ij as well as from the discreteness and the 
finite number of individuals. 

To investigate the behavior of the interacting particle 
system and to compare it with the predictions of the 



straightforward algorithm, at each Monte Carlo (MC) step, 
one determines (via a random number) which reaction (ex- 
change, selection, or reproduction) occurs next. Reproduc- 
tion happens with probability /i/(/i + a + e), selection with 
a/di + a + e), and exchange events occur with probability 
e/{fi + (7 + e). Then, a random pair of nearest neighbors is 
selected and the chosen type of interaction (reproduction, 
selection or exchange) is performed, if the move is allowed. 
In our situation, a more efficient algorithm inspired by |Gille-"] 
spie ( 1976 1977 1 can be implemented. Namely, in the con- 
tinuum limit N —^ oo, the exchange rate e becomes large 
compared to the selection and reproduction rates, fj, and a. 
Thus, a large number of exchange events occurs between 
two reactions ([iJ. Indeed, the probability P of having E 
exchanges between two subsequent May-Leonard reactions 
reads 

E / 



PiE) 



(3) 



Hereby, the first factor on the right hand side denotes the 
probability of subsequently drawing E exchange events, and 
the second factor is the probability that the next event is 
either a selection or reproduction process. To efficiently 
take this high number of exchanges occurring between selec- 
tion/reproduction processes into account, at each MC step, 
we draw the number of such exchange events from the prob- 
ability distribution ^ . This number of random exchanges 
is performed, and then followed by one of the May-Leonard 
reactions Q. 

All results we present from lattice simulations were ob- 
tained starting from a random initial distribution of indi- 
viduals and vacancies: the probability for each site to be 
in one of the four possible states (i.e. A,B,C or 0) has 
been chosen to coincide with the value of the (unstable) 
internal fixed point of the rate equations ([4]). Thereafter, 
without loss of generality (see Section 3.6), we often con- 
sider equal selection and reproduction rates, which we set 
to one (thereby defining the time scale), i.e. fi = a = I. In 
this case, all four states initially occur with equal probabil- 
ity 1/4. Generally, after a short transient, a reactive steady 
state emerges. Hereby, the discussion of the stability of the 
reactive steady states should be dealt with some care be- 
cause the only absorbing states are those where the system 
is uniformly covered by only one subpopulation. Reactive 
states are not stable in a strict sense because they can be 
spoilt by chance fiuctuations (the system is large but fi- 
nite) which drive the system into one of the absorbing uni- 
form states. However, the typical (average) waiting time T 
to reach the extinction of two subpopulations is extremely 
long for large systems and it diverges with increasing system 



SPDE (30), we have carried out stochastic simulations of size (Reichenbach et al. 2007a I, implying the existence of 



the model on a square lattice with periodic boundary con- 
ditions and of size ranging from L = 50 up to L = 1000 
sites. In the following, we always consider the system in 
two spatial dimensions, d = 2. At each simulation step, a 
randomly chosen individual interacts with one of its four 
nearest neighbors, being also randomly determined. In a 



super-persistent transients (Hastings 20041. To rationalize 



this point, here we follow (Reichenbach et al. 2007a I where 



we have proposed to characterize the stability of the reac- 
tive steady states by comparison with the average extinction 
time obtained for the marginally stable version of the sys- 
tem, where T oc N (Reichenbach et al. 20061. Therefore, 
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when the extinction time grows faster than hnearly with the 
system size, the reactive steady state is deemed to be stable 
on ecologically reasonable time scales. 

In the continuum limit, the system is described in terms 



L = 100 



L = 200 



ociV- 



of the SPDE ( 30 ) , characterized by white noise of strength 
'^/"^ (see Section 3.5), which are derived from the mas- 



ter equation via a system-size expansion as explained in 
Appendix D. To compare this approach with the results of 



the lattice system, we have numerically solved Eqs. (30 1 us- 



ing open software from the XMDS project (Collecutt and 



Drummond 2001 XmdS I. For specificity, we have used spa- 



tial meshes of 200 x 200 to 500 x 500 points, and 10,000 steps 
in the time-direction. Initially, we started with homoge- 
neous densities, corresponding to the values of the unstable 
internal fixed point of the RE Q, as in the lattice simu- 
lations. It is worthwhile noticing that noise terms of the 
SPDE have to be treated with special care as they may oc- 
casionally lead to the nonphysical situation of negative den- 
sities (or densities exceeding the maximal value of 1) if the 
system is close to the phase space boundaries. This prob- 
lem is most pronounced when additive noise is used, while 
in our case, noise is multiplicative and vanishes near the 
boundaries. Still, due to discretization effects, nonphysical 
situations of negative densities may arise; in our simulations, 
we have discarded such rare events. 

2.1 Spiral structures in the continuum 
limit 

To study the system's behavior in the approach of the con- 
tinuum limit, we have kept the diffusion fixed at a value 
D = \ X 10^^, and systematically varied the system size N 
(and henceforth the exchange rate e, which reads e = 2DN 
in two spatial dimensions). In Fig. [3] we report typical long- 
time snapshots of the reactive steady states for various val- 
ues of the exchange rate and different system sizes. In small 
lattices, e.g. L = 100, we observe that all subpopulations 
coexist and form clustering patterns of characteristic size. 
The spatio-temporal properties of the latter do not admit 
simple description and appear to be essentially dominated 
by stochastic fluctuations. When the size of the lattice is 
large, e.g. L = 300 — 500, the three populations also co- 
exist in a reactive steady state. However, in stark contrast 
from the above scenario, we now find that individuals self- 
organize in an entanglement of rotating spiral waves. The 
properties of these spirals, such as their wavelength, their 
frequency and the spreading velocity, are remarkably regu- 
lar and will be studied below. 

The asymptotic approach towards the continuum limit, as 
illustrated in the snapshots of Fig. [3] can be rationalized by 
considering the typical size of the patterns. The latter are 
obtained from the computation of spatial correlation func- 
tions (see next Subsection for the definition of correlation 
functions) and the ensuing correlation length ^com which is 
the length at which the correlations decay by a factor 1/e 
from their maximal value). This quantity gives the aver- 
age spatial extension of the patterns. In Fig. [4] we report 




Figure 3: Approach of the continuum limit. We show snap- 
shots of the reactive steady state of the stochastic system, for 
D = 1 X 10~^, fi — a = 1, and different system sizes. Each color 
(level of gray) represents a different species (black dots denote 
empty spots). The lattice sizes are L = 100 (e — 0.2) in the 
top left panel, L = 200 (e = 0.8, top right), L = 300 (e = 1.8, 
bottom left), and L = 500 (e — 5, bottom right). Increasing the 
system size {D is kept fixed), the continuum limit is approached. 
Random patterns appear for small systems (L = 100), while en- 
tangled spiral waves emerge when L is raised and clearly emerge 
in large systems {L — 500). 



0.25 




Figure 4: The typical size of the patterns, as described by the 
correlation length £corr. We show its dependence on the exchange 
rate e, for fixed D — 5 x 10~^, fi = a = 1, and different system 
sizes A''. For large A'', i.e. large e, we expect the system to be 
well described by SPDE, the correlation length of the latter is 
depicted as a dashed line. Surprisingly, the correlation length 
already agrees excellently with this continuum model for e > 5. 
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D = 1 X 10 




Figure 5: The reactive steady states. We show snapshots emerging in simulations of the interacting particle system |T| (top row) 
and obtained by solving the SPDE (301 (bottom row). Each color (level of gray) represents a different species (black dots denote 
empty spots). From left to right, the diffusion constant is increased from D = Ix 10~® to D = 3 x 10"'*. The latter value is slightly 
below the critical threshold above which the spiral structures can no longer fit within the system ( [Reichenbach et al.[ pOOTa] ); see 
text. The system sizes used in the stochastic simulations are L = 1000 in the left two panels, L — 300 for that at right, and L = 500 
for the other two (middle). The selection and reproduction rates are chosen as a = /i = 1. 



£corr as obtained for systems of different sizes, i.e. for var- 
ious values of e. For small systems (low exchange rate e), 
the pattern size is considerably larger than in the contin- 
uum limit (dashed line obtained from SPDE, see below). 
Increasing the system size (larger values of e), the contin- 
uum limit is approached. Remarkably, there is a striking 
agreement between the results from the lattice simulations 
and the SPDE already for e > 5. Hereby, ii — a — 1, 
such that the continuum limit is already closely approached 
when e is of the same order as the rates for the selection 
and reproduction events, yet larger. This result is also ap- 
parent in Fig. [Sj where for L = 500 and e — 5, the system 
already exhibits regular spiral waves. It follows from this 
discussion that the results obtained in the continuum limit 
(derived assuming e oo), actually have a broader range 
of validity and allow to aptly describe the interacting parti- 
cle system already for e finite and of the same order of (yet 
larger than) fi and a. A comparable influence of short-range 
mixing has also been reported recently for a predator-prey 
- or host-pathogen - model, where a short-range exchange 
process with finite rate has been shown to crucially affect the 
fate of the system (absorbing or coexistence state) through 
a (first-order) phase transition ( Mobilia et al.[ [2006 1 . Fur- 
thermore, the smooth domain boundaries caused by mixing 
and the emerging spiral waves are similar to the spatial pat- 



terns investigated in (Szabo and Szolnoki 2002) that arise 



from slow cyclic dynamics combined with Potts energy. The 
authors of this study have analyzed the resulting spirals by 
considering the vortex density and average length of vortex 
edges. 

In Fig. [5] we report snapshots of the reactive steady state 



obtained from the stochastic simulations (left column) and 
as predicted by the SPDE (right column) for different values 
of the diffusion constant D. In Fig. [5] panels in the same 
row have been obtained for the same parameters (for lat- 
tices of size L ~ 300,500 and 1000). We observe an excel- 
lent qualitative and quantitative agreement between both 
descriptions, which yield the formation of rotating spiral 
waves whose typical sizes and wavelengths manifestly co- 
incide (see below). When the diffusion constant D is in- 
creased, the size of the spirals increases, too. With help of 



the underlying SPDE ( 30 1 , this observation can be ratio- 



nalized by noting that the size of the spatial structures is 
proportional to ^/D. This scaling relation stems from the 



fact that spatial degrees of freedom only enter the Eqs. ( 30 ) 



through the diffusion term DA, where the Laplacian oper 
ator A involves second-order spatial derivatives. Therefore, 
rescaling the spatial coordinates by a factor l/VlD makes 
the diffusive contribution independent of D, with the en- 
suing scaling relation. As we have shown in (Reichenbach 



et al. 



2007a), the three subpopulations can coexist only up 
to a critical value of the diffusion rate. Above that thresh- 
old, the spirals outgrow the system and there is extinction 
of two populations: only one subpopulation (at random) 
survives and covers the entire lattice. 



2.2 Correlations 

The comparison of snapshots obtained from lattice simula- 
tions with the numerical solutions of the SPDE reveals a 
remarkable coincidence of both approaches (see Fig. [5| . Of 
course, due to the inherent stochastic nature of the inter- 
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Figure 6: Correlation functions. The spatial correlation gAA{r) 
as function of r in the reactive steady state is shown. We report 
results obtained from stochastic simulations (red circles, for a 
lattice of linear size L — 1000) and numerical solutions of the 
SPDE (30 1, blue squares, and notice an excellent agreement. In 
both cases, results have been obtained for a value D = 3 x 10~® 
and ^ = a — 1. The typical correlation length £corr as a function 
of the diffusion constant D is shown in the inset (on a double 
logarithmic scale). The scaling relation Icon ~ VT^, indicated 
by a black line, is clearly confirmed. We have also reported the 
results for the static correlation function gAA{r) of the patterns 
predicted by the deterministic PDE (green triangles); see text. 
The latter are found to be markedly less damped than those 
arising in the stochastic descriptions of the system. 



acting particle system, the snapshots do not match exactly 
for each realization. To reach a quantitative assessment 



on the validity of the SPDE (30 1 to describe the spatio- 
temporal properties of the system in the continuum limit, 
we have computed various correlation functions for the sys- 
tem's steady state. The attainment of the steady state is as- 
sessed by computing the long time evolution of the densities 
and various snapshots as those of Fig. [5] When the densities 
are found not to fluctuate significantly around their average 
values and the snapshots display statistically the same ro- 
bust features at various times (typically t ^ 100—1000), the 
system is considered to be settled in its (reactive) steady 
state. 

We first consider equal-time correlation functions, which 
yield information about the size of the emerging spirals. 
As an example, we focus on the correlation gAAil^ — r'\) 
at r and r' of the subpopulation A, gAA{\f — r'\) — 
{a{r,t)a{r' ,t)) — {a{r , t)) {a{r' , t)) . Here, the brackets (...) 
stand for an average over all histories. In the steady state, 
the time dependence drops out and, because of translational 
and rotational invariance, the latter depends only on the 
separating distance |r — r'|. In Fig. [6] we report results for 
gAA obtained from lattice simulations (red circles) and from 
numerical solutions of the SPDE ( pO] ) (blue squares) , finding 
an excellent agreement between them. When the separat- 
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Figure 7: Autocorrelation. We show the correlation gAA{t) as a 
function of time t. Results from stochastic simulations (red/gray 
line) are compared with those obtained from the numerical solu- 
tions of the SPDE (blue squares), as well as with those computed 
from the deterministic PDE (green triangles). All results are in 
excellent agreement with each other and are characterized by 
oscillations at frequency fi""™ ^ 0.103 (for fi = a = 1); the 
latter is independent from the value of the diffusion D. These 
oscillations reflect the rotation of the spiral waves. The results 
from the SPDE and deterministic PDE have been obtained using 
D = 10~^, while stochastic simulations have been performed on 
a lattice of length L = 300 with D = 10"''. 



ing distance vanishes, the correlation reaches its maximal 
value and then decreases, exhibiting (damped) spatial os- 
cillations. The latter reflect the underlying spiralling spa- 
tial structures, where the three subpopulations alternate in 
turn. Damping results from the averaging over many small 
spirals. As described in the previous subsection, the correla- 
tion functions are characterized by their correlation length 
^corr, which conveys information on the typical size of the 
spirals. In the inset of Fig. [6] we show the dependence of the 
correlation length on the diffusion rate D in a double loga- 
rithmic plot which confirms the scaling relation ^corr V~D, 
also inferred from general considerations. 

We now consider the time dependence of the correlation 
functions and study the autocorrelation gAAi\t — t'\) of sub- 
population A at times t and t' , for a fixed spatial position. 
This quantity is given by gAA{\t — t'\) — {a{r,t)a{r,t')) — 
{a{r,t)){a{r,t')) and only depends on the time difference 
\t — t'\. Both lattice simulations and SPDE ( [SO] ) yield os- 
cillating correlation functions, as shown in Fig. I7j This pe- 
riodic behavior, with a frequency numerically found to be 
j-^num ^ 0.103 (for (T — /i = 1), stcms from the rotational 
nature of the spiral waves and is independent of the diffusion 
constant D. Below, this value is compared with an analyt- 
ical prediction inferred from a deterministic description of 
the spatial system. In the time intervals which we have in- 
vestigated, t ~ 1,000, the oscillations, as reported in Fig. [7] 
are undamped. Therefore, on this time-scale, the position of 
the spirals' vortices is stable in the steady state and not in- 
fluenced by noise. On larger time-scales, however, we expect 



the vortices to perform random walks (see (Cross and Ho- 



henberg 19931 for a general discussion as well as (Tainaka 
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Figure 8: The role of noise. Left: numerical solution of the 
SPDE | |30[ ), starting from an homogeneous initial state (all sub- 
populations have an initial density 1/4). Each color (level of 
gray) represents a different species (black dots denote empty 
spots). Right: if we ignore the noise terms of the SPDE 1301) 



and study the resulting deterministic PDE ( 23 1, the steady state 



depends on the initial configuration. Here, this is illustrated by 
a snapshot of the steady state evolving from the initial condi- 
tion: a(0) = a* + cos(27rxy)/100, 6(0) — b* , c(0) = c*, see 
text. In both panels, the diffusion constant is _D = 3 x 10~^ and 
cr = = 1. 



1994 Nishiuchi et al. 2008') for investigations of vortex dy- 



namics in rock-paper-scissors models), with associated vor- 
tex annihilation and creation processes. Studies exploring 
such a behavior are promising for further broadening the un- 
derstanding of stochastic effects on nonequilibrium steady 
state. In the following, we consider another footprint of 
stochastic fluctuations: we show how the latter significantly 
influence the steady state of the system when starting from 
homogeneous initial densities. 

2.3 The role of stochasticity 

So far, we have considered stochastic descriptions, relying 
on lattice simulations and the stochastic partial differential 
equations (30). In the latter, the noise terms are propor- 
tional to i.e. vanishing in the limit of iV — > oo. 
It is therefore legitimate to ask whether it is possible to 
simply neglect noise and describe the system in terms of de- 
terministic partial differential equations (PDE) ( pS] ) given 
in Subsection |3.4[ To address this question and to reach 
a better understanding of the effects of internal noise, we 



have numerically solved the deterministic PDE ( 23 ) for var- 
ious initial conditions. Of course, the latter have always 
to be spatially inhomogeneous, otherwise spatial patterns 
cannot emerge. Starting from a spatially inhomogeneous 
distribution of the populations, the deterministic equations 
are found to evolve towards a reactive steady state also char- 
acterized by the emergence of spirals. We have numerically 
checked that these spirals' wavelengths A are the same as 
those obtained in the stochastic lattice simulations and for 



the solutions of the SPDE (30) 



However, there are major differences between the deter- 
ministic and the stochastic descriptions of the spatially ex- 
tended system. Concerning the SPDE, when the initial den- 



sities are the ones corresponding to the unstable internal 
fixed point of the rate equations ([4| , even without initial in- 
homogeneities, an entanglement of spiral waves arises in the 
course of time evolution only due to noise. We have numer- 
ically found that the latter is characterized by a universal 
vortex density of about 0.5 per square wave length. For the 
deterministic PDE, spatially homogeneous initial conditions 
do not yield spiralling patterns. When starting with initial 
spatial inhomogeneities, the density of the latter sensitively 
determines the density of spirals, which can be much lower 
than in the stochastic situation. As an illustration, in Fig. [8] 
we compare snapshots of the steady state of the SPDE and 
of the PDE. For the latter, we have chosen the initial con- 
dition a{t = 0) = a* +cos(27ra;y)/100, &(0) = &*, c(0) = c*, 
just adding a small local perturbation to the value of the 
homogeneous fixed point (a*,&*,c*). While a large number 
of spirals cover the system in the stochastic case (Fig. [s] 
left), only four spirals appear in the deterministic situation; 
Fig. [8] right. This differe nee is also manifest when one con- 
siders the spatial dependence of the correlation functions, 
as shown in Fig. [6] These quantities share the same max- 
ima and minima for the stochastic and deterministic de- 
scriptions, which can be traced back to the fact that spirals 
have the same wavelengths, not affected by the noise. How- 
ever, the correlations obtained from lattice simulations and 
from the SPDE are much more strongly damped than in the 
case of a deterministic PDE. This stems from the fact that 
stochasticity effectively acts as an internal source of spa- 
tial inhomogeneities (randomly distributed), resulting in a 
larger number of spirals (of small size). The agreement be- 
tween the temporal dependence of the correlation functions 
in the deterministic and stochastic descriptions (see Fig. [7| 
is another consequence of the fact that the spirals' frequency 
is not affected by the noise. 



We can now ask what happens if one solves the SPDE ( 30 ) 



with inhomogeneous initial conditions. To answer this ques- 
tion, we have systematically studied the emerging steady 
state upon varying the strength of noise and the degree of 
spatial inhomogeneity of the initial configuration. When 
these effects are of comparable intensity (i.e. when the noise 
strength is of the same order as the initial deviations of the 
densities from spatial homogeneity) , the steady state is still 
driven by noise and gives rise to an entanglement of small 
spirals. In this situation, the universal density of 0.5 spiral 
vortices per square wavelength arises. On the other hand, 
if the degree of inhomogeneity of the initial state is signifi- 
cantly higher than the noise level, the former dominates the 
systems' behavior and the density of spirals is determined 
by the spatial structure of the initial configuration, as in 
Fig. [s] (right). Therefore, if no initial inhomogeneities are 
present, stochasticity acts as a random source of inhomo- 
geneities leading to the emergence of an entanglement of 
stable spiral vortices. The latter are stable against stochas- 
tic effects, as reflected by the undamped temporal oscilla- 
tions of the autocorrelation function of Fig. [7] On the other 
hand, if initial inhomogeneities exceed the noise level, they 
are responsible for the formation of vortices before noise can 
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influence the system. 

From the above discussion, we infer that the range of ini- 
tial conditions where noise influences the system's fate and 
leads to the universal density of about 0.5 spiral vortices 
per square wavelength is rather small. It contains initial 
densities that lie around the values of the unstable reactive 
fixed point of the rate equations ^ , with spatial variations 
whose amplitude does not exceed the noise level. Other 
initial conditions, with more strongly pronounced spatial 
inhomogeneities, lead to spiral vortices that are determined 
by the initial inhomogeneities, and not stochastic effects. 
The latter behavior is remarkable, as it corresponds to a 
memory of the system: The information about the initial 
state, i.e. its spatial inhomogeneities, is preserved in the 
reactive steady state where it manifests itself in the posi- 
tion of the spirals' vortices. In the time interval that we 
have considered (t ~ 1000), the latter are stable during the 
temporal evolution and noise does not affect their location. 
This feature is intimately connected to the nonequilibrium 
nature of the system. Indeed, in noisy equilibrium systems, 
the steady state does not convey information about the ini- 
tial conditions. This stems from the standard assumption 
of ergodicity, as first formulated in the works of Boltzmann 
and Gibbs fBoltzmannI 119641 iGibbsl 119021) . Only a genuine 



nonequilibrium steady state can display memory of the ini- 
tial state. 

Finally, let us comment on the nature of the noise en- 
countered in our analysis. In the SPDE (30 1, the noise is 



multiplicative because its strength depends on the densi- 
ties of the subpopulations. Approximating the latter by 
additive white noise [e.g., substituting (a, 6, c) by the values 
(a* , 6* , c* ) in the expressions of the noise contributions] , we 
essentially obtain the same results as with multiplicative 
noise. We understand this behavior as stemming from the 
fact that the main infiuence of noise is to spatially and ran- 
domly perturb the system's initial state. Hence, as long its 
intensity scales like N~^^'^ (weak noise), the fact that noise 
is either multiplicative or additive has no significant impact 
on the system's fate. 

2.4 The spirals' velocities, wavelengths, 
and frequencies 

Above, we have found that characteristic properties of the 
emerging spiral waves, like their wavelength and frequency, 
are unaffected by noise. To compute these quantities an- 
alytically, it is therefore not necessary to take noise into 
account, and we may focus on the study of the determinis- 
tic PDE (23 1. In Subsection [3^ we show how the dynamics 



of the latter is essentially captured by an appropriate com- 



plex Ginzburg-Landau equation (CGLE) , given by Eq. ( 32 1 
for the case under consideration here. The CGLE (32 1 al- 



lows to derive analytical results for the emergence of spiral 
waves, their stability and their spreading velocity, as well as 
their wavelength and frequency. We detail these findings in 
Subsections |3.7| and |3.8| Here, we assess the accuracy and 
validity of these analytical predictions by comparing them 
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Figure 9: Spreading velocity. We report the dependence of 
front velocity v* (rescaled by a factor VTj) on the reproduction 
rate fi. The time scale is set by keeping a — 1. In red (full 
line), we report the analytical predictions (371 obtained from 



the CGLE, which are compared with numerical results (black 
dots). The latter are obtained from the numerical solutions of 



the SPDE (301 



with values obtained from the numerical solutions of the 



SPDE (301 



Let us first consider the spreading velocity v* of the 
emerging wave fronts. The analytical value, inferred from 



the CGLE and derived in Subsection |3J| [see Eq. p7|)] 
reads v* = 2^/ciD, where ci = fia/[2{3fj, + a)] is a coef- 
ficient appearing in the CGLE ( 32 1 . In numerical compu- 



tations, the front velocity is obtained from the wavelength 
A and the frequency f2 of the emerging spirals. Namely, 
the wavelength A""™ can be inferred from snapshots (as in 
Fig. [5| , and the frequency fi""™ is computed from the os- 
cillations of the autocorrelation (as in Fig.]?]). The velocity 
then follows via w"""^ = A""'"r2"^'"/27r. As the wavelength 
is proportional to y/D and the frequency does not depend 
on the diffusion constant, one can easily check that the re- 
lation w""'" = A™'"r2"""V27r confirms that u"^" ~ /D, as 
m Eq. ([37]). In Fig. [9] we compare the analytical prediction 
(37) for V* with results obtained from the numerical solu- 
tion of the SPDE ( [30| , as function of the reproduction rate 
/i (setting cr = 1, we fix the time scale), and find a good 
agreement. On the one hand, for small values of fi (much 
lower than the selection rate, /i ^ 1), reproduction is the 
dominant limiter of the spatio-temporal evolution. In the 
limit /Lt — > 0, the front velocity therefore only depends on 
/z. From dimensional analysis, it follows v* ~ as also 
confirmed by the analytical solution Eq. ( 37 1 . On the other 
hand, if reproduction is much faster than selection, /i 3> 1, 
the latter limits the dynamics, and we recover v* ^ \Ja. 
In Fig. [9] as (7 = 1, this behavior translates into v* being 
independent of in this limit. While the numerical and 
analytical results coincide remarkably for low reproduction 
rates (i.e. /i < 0.3), systematic deviations (« 10%) appear 
at higher values. As an example, when selection and re- 
production rates are equal, cr = /x = 1 (as was considered 
throughout the last section), we have numerically found a 
velocity t;™™ « 0.63v^, while Eq. (37) yields the analyti- 
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Figure 10: The spirals' wavelength. We show the functional 
dependence of the wavelength A on the rate \i (with a = 1), 
and compare numerical results (black circles) , obtained from the 



numerical solutions of the SPDE (30 1, to analytical predictions 
(red line). The latter stem from the CGLE and are given by 
Eq. (42 1. They differ from the numerics by a factor of 1.6, see 
text. Adjusting this factor, c.f. the blue line, the functional 

dependence is seen to agree very well with numerical results. 3.1 Rate equations 



approaching the boundaries of the phase space, while the 
dynamics underlying the CGLE is characterized by limit 
cycles (usually distant from the edges of the phase space) 
resulting from a (supercritical) Hopf bifurcation. 

3 Theory 

In the following, we present and discuss an analytical 
approach built on stochastic partial differential equations 
(SPDE) arisen in the proper continuum limit of large sys- 
tems with mobile individuals. Applying the theory of in- 
variant manifolds and normal forms to the nonlinear parts 
of these equations, we show that the SPDE fall into the 
universality class of a complex Ginzburg-Landau equation 
(CGLE) . We derive the latter and employ it for a quantita- 
tive determination of properties of the spiral waves. 



cal result v* = y/Dj2 « 0.71/D. 

Concerning the spirals' wavelengths and frequencies, in 
Subsection |3.8[ we analytically infer predictions from the 

We have checked 
In Fig 



CGLE (321 given by Eqs. (41 1 and (42 1 



The deterministic rate equations (RE) describe the tempo- 
ral evolution of the stochastic lattice system, defined by the 
reactions ([T]), in a mean- field manner, i.e. they neglect all 
spatial correlations. They may be seen as a determinis- 
tic description (for example emerging in the limit of large 
system sizes) of systems without spatial structure, such as 



these results against numerical computations 
the analytical estimates for the wavelength A are compared 
with those obtained from the numerical solution of the 



10] Moran processes (Moran 1958 Traulsen et al. 2005 2006) 



or urn models (Feller 1968 Reichenbach et al. 20061. The 



SPDE (30 1 for different values of the reproduction rate /i. 



We notice that there is an excellent agreement between ana- 
lytical and numerical results for the functional dependence 
of A on /i. For low reproduction rates (/i <C 1) we have 
A ^ 1 1 y/Ji, while when reproduction occurs much faster 
than selection (/x 3> 1), the dynamics is independent of /i 
and A ~ l/v^- We have also found that the analytical 
result predicts an amplitude of A which exceeds that ob- 
tained from numerical computations by a constant factor 
K, 1.6, taken into account in Fig. [TOl W e attribute this de- 
viation to the fact that the CGLEp2|) (stemming from the 



study of the RE is the ground on which the analysis of the 
stochastic spatial system is built. In particular, the proper- 
ties of the RE are extremely useful for the derivation of the 



system's SPDE (jSOf and CGLE ([321. 

Let a, 6, c denote the densities of subpopulations A, B, 
and C, respectively. The overall density p then reads p = 
a + b ~\- c. As every lattice site is at most occupied by one 
individual, the overall density (as well as densities of each 
subpopulation) varies between and l,i.e. 0<p<l. With 
these notations, the RE for the reaction ([!]) are given by 



normal form (13l) describes a dynamics exhibiting a limit 
cycle, while the full May-Leonard rate equations ^ are 
characterized by heteroclinic orbits. The correct functional 
dependence of the wavelength A on the reproduction rate 
p is therefore especially remarkable. Elsewhere it will be 
shown that in the presence of mutations, inducing a limit- 
cycle behavior, the description of the emerging spiral waves 



dta 
dtb 
dtc 



a[p{l ~ p) 
h[p{l-p) 
c[p{l - p) 



ac\ 



ah] . 



(4) 



They imply the following temporal evolution of the total 
density: 



dtp — /ip(l — p) ~ a{ab + hc + ac) . 



(5) 



in terms of CGLE (32 1 becomes fully accurate. 



For the spirals' frequency, we analytically obtain Q. — 
n{q'°^) ^ LU + 2(ci/c3) (l - yiTcf), see Subsection 
As already inferred from numerical simulations (Sec. Ill 



These equations have been introduced and investigated 



3.8 



by May and Leonard (19751. In the following, we review 



f2 does not depend on the diffusion D. Quantitatively, and 
as an example for p = a = I, we obtain the analytical 
prediction « 0.14, which differs by a factor « 1.4 from 
the numerical value fi""™ « 0.103 found in Fig. [7] As for 

the wavelength, this difference stems from the fact that the subpopulation, {al,b2,C2) 
May-Leonard rate equations ^ predict heteroclinic orbits 



some of their properties. 

Equations (|4| possess four absorbing fixed points. One 
of these (unstable) is associated with the extinction of all 
subpopulations, {al,bl,cl) — (0,0,0). The others are het- 
eroclinic points (i.e. saddle points underlying the hetero- 
clinic orbits) and correspond to the survival of only one 



(1,0,0), (a^,&*,c;) = (0,1,0) 
and (a|,6|,c|) = (1,0,0), shown in blue (dark gray) in 
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Fig. [2] In addition, there exists a fixed point, indicated 
in red (gray) in Fig. |2] where all three subpopulations coex- 
ist (at equal densities), namely {a*,b*,c *) = ^^^^ (1,1,1). 
For a non- vanishing selection rate, a > 0, |May and Leonard] 
( 1975 1 showed that the reactive fixed point is unstable, and 



the system asymptotically approaches the boundary of the 
phase space (given by the planes a = 0, & = 0, and c = 0). 
There, they observed heteroclinic orbits: the system oscil- 
lates between states where nearly only one subpopulation 
is present, with rapidly increasing cycle duration. While 
mathematically fascinating, this behavior was recognized to 



be unrealistic (May and Leonard 19751. For instance, in a 



biological setting, the system will, due to finite-size fluctua- 
tions, always reach one of the absorbing fixed points in the 
vicinity of the heteroclinic orbit, and then only one popula- 
tion survives. 

Linearization of the RE Q around the reactive fixed 
point leads to dtx — Ax with the vector x — {a — a* ,b — 
b* ,c— c*)"^ and the Jacobian matrix 



A 




(6) 



As this matrix is circulant, its eigenvalues can be obtained 
from a particularly simple general formula (see e.g. (Hof- 



bauer and Sigmund[ |1998) ); they read: 



Ao 
Ai 

A, 



1 /icr 



2 3fi + a 

1 /XtT 

2 3n + a 



[1 -I- Vii] 
[1 - Vsi] 



(7) 



This shows that the reactive fixed point is stable along the 
eigendirection of the first eigenvalue Aq. As elaborated be- 



low, there exists an invariant manifold (Wiggins 1990) (in- 



cluding the reactive fixed point), that the system quickly 
approaches. To first order such a manifold is the plane nor- 
mal to the eigendirection of Aq. On this invariant manifold, 
flows spiral away from the reactive fixed point, which is an 
unstable focus, as sketched in Fig. |2] (blue trajectory). 

The linear stability analysis only reveals the local stability 
of the fixed points. The global instability of the reactive 
fixed point is proven by the existence of a Lyapunov function 



£ ( jHofbauer and Sigmund[ |1998[ [May and Leonard] |1975[ ): 



C = 



ibc 



(8) 



In fact, using Eqs. ^ and ([5]), the time derivative of L is 
found to be always non-positive, 

dtC^ -]^ap-^abc[{a~bf + {b~cf + {c-af] < . (9) 

We note that dtC vanishes only at the boundaries (a — 
0, 6 = or c = 0) and along the line of equal densities, 
a — b—c. The latter coincides with the eigendirection of Aq, 
along which the system approaches the reactive fixed point. 



However, on the invariant manifold we recover dtC < 0, 
corresponding to a globally unstable reactive fixed point, as 
exemplified by the trajectory shown in Fig. |2] 

As the RE ([4| have one real eigenvalue smaller than zero 
and a pair of complex conjugate eigenvalues, they fall into 
the class of the Poincare-Andronov-Hopf bifurcation, well 
known in the mathematical literature (Wiggi ns, ^1990) . The 
theory of invariant and center manifolds allows us to recast 
these equations into a normal form. The latter will turn 
out to be extremely useful in the derivation of the CGLE. 
In the following, we derive the invariant manifold to second 
order as well as the normal form of the RE. 

3.2 Invariant manifold 

An invariant manifold is a subspace, embedded in the phase 
space, which is left invariant by the RE (|4|, i.e. by the de- 
terministic dynamics. In the phase space, this means that 
fiows starting on an invariant manifold always lie and evolve 
on it. Here, we consider a two-dimensional invariant man- 
ifold associated with the reactive fixed point of the RE ([4| 
onto which all trajectories (initially away from the invariant 
manifold) decay exponentially quickly (see below). Upon 
restricting the dynamics to that appropriate invariant man- 
ifold, the system's degrees of freedom are reduced from three 
to two, which greatly simplifies the mathematical analysis. 

To determine this invariant manifold, we notice that the 
eigenvector of the eigenvalue Aq < at the reactive fixed 
point is a stable (attractive) direction. Therefore, to lowest 
order around the reactive fixed point, the invariant man- 
ifold is simply the plane normal to the eigendirection of 
Aq. As well known in the theory of dynamical systems, 
beyond first order, nonlinearities affect the invariant man- 
ifold, which is therefore only tangent to this plane. To in- 
clude nonlinearities, it is useful to introduce suitable coor- 
dinates (y^, ?/B, j/c) originating in the reactive fixed point. 
We choose the yc-axis to coincide with the eigenvector of 
Aq, such that the plane yc = is the invariant manifold 
to linear order, onto which the fiows of ^ relax quickly. 
The coordinates ua and i/b are chosen to span the plane 
normal to the axis ycj forming an orthogonal set. Such co- 
ordinates y = {yA,yB,yc) are, e.g., obtained by the linear 
transformation y = Sx, with the matrix S given by 



S 




(10) 



where x — {a — a*,b — b*,c — c*)^ denotes the deviation 
of the densities from the reactive fixed point. The coordi- 
nates {yA^yBTVc) are shown in Fig. [2] and their equations 
of motion are given in Appendix A. 

The May-Leonard model is symmetric under cyclic per- 
mutations of A, B, and C. In the y-coordinates, each of 
these permutations translates into a rotation of 27r/3 around 
the yc axis. The equations of motion refiect this symmetry 
of the system, as can be checked explicitly in the Eqs. (44) 
in Appendix A. 
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To parameterize the invariant manifold sketched in Fig. [2] 
we seek a function G{yA,yB), with yc = G{yA,yB)- If all 
nonlinearities of the RE are taken into account, this is a 
very complicated problem. However, for our purpose it is 
sufficient to expand G to second order in yA,yB- As the 
invariant manifold is left invariant by the RE, by definition, 
G must obey 



the new variables: 



dG 

dtG{yA{t),yB{t)) = ^ — dtyA 
oyA 



dG 
dVB 



dtyB = dtyc 



(11) 

To linear order in yA and yB, '^'^ simply have G = and re- 
cover yc — 0, corresponding to the plane normal to the 
j/c-direction. We have anticipated this result above: to 
first order, the invariant manifold coincides with this plane, 
and is tangential to it when higher orders are included. To 
second order, only linear terms of dtyA,dtyB contribute to 
Eq. ([TT]). The latter are invariant under rotations in the 
(j/yi, yB)-plane, and G must obey the same symmetry. It is 
therefore proportional to y\ + y^. After some calculations, 
detailed in Appendix B, one obtains: 



yc = G{yA,yB) 



a 3^ - 



4/i 3/i + 2f7 



{y\ + yl) + o{y^). (12) 



The comparison of this expression for the invariant mani- 
fold, valid to second order, with the numerical solutions of 
the RE Q (which should, up to an initial transient, lie on 



the invariant manifold) confirms that (12 1 is an accurate 



approximation, with only minor deviations occurring near 
the boundaries of the phase space. 

3.3 Normal form 

Nonlinear systems are notably characterized by the bifur- 
cations that they exhibit (Wiggins 1990). Normal forms 



are defined as the simplest differential equations that cap- 
ture the essential features of a system near a bifurcation 
point, and therefore provide insight into the system's uni- 
versal behavior. Here, we derive the normal form associ- 
ated with the RE ^ of the May-Leonard model and show 
that they belong to the universality class of the Hopf bifur- 
cation (Wiggins 19901. Below, we demonstrate that this 



property allows to describe the system in terms of a well- 
defined complex Ginzburg-Landau equation. 

Restricting the (deterministic) dynamics onto the invari- 
ant manifold, given by Eq. (12 1, the system's behavior can 

Here, we choose to 
with the resulting 



be analyzed in terms of two variables, 
express yc as a function of yA and yB, 
rate equations (up to cubic oder) given in Appendix A. The 



latter can be cast into a normal form (see (Wiggins 19901 



Chapter 2.2) by performing a nonlinear variable transfor- 
mation y z which eliminates the quadratic terms and 
preserves the linear ones (i.e. y and z coincide to linear or- 
der). As an ansatz for such a transformation, we choose the 
most general quadratic expression in y for the new variable 
z. Details of the calculation can be found in Appendix C. 
Here we quote the result for the normal form of the RE in 



dtZA 
dtZB 



CiZA - 
ClZB 



- UJZA 



C2 {ZA + CgZs) {z\ + z]j) + o{z^) , 

C2{zB - C3Za){za + z%) + o{z^) . 



In these equations, 



-\/3 /i(T 



(13) 



(14) 



is the (linear) frequency of oscillations around the reactive 
fixed point. The constant 



1 



Cl 



/icr 



2 3/i -f cr ' 



(15) 



gives the intensity of the linear drift away from the fixed 
point, while 



C2 



C3 



cr(3/Lt-h cr)(48/x-f llcr) 



(16) 
(17) 



are the coefficients of the cubic corrections. 

To gain some insight into the dynamics in the normal 



form, it is useful to rewrite ( 13 1 in polar coordinates (r, i 
where za 



r cos (p, zb = rsiTL (, 

dtr = r[ci 
dtO = —u! - 



. This leads to 



(18) 



These equations only have a radial dependence, which 
clearly reveals a polar symmetry. They predict the emer- 
gence of a limit cycle of radius r — \fc\fc2 and therefore 
fall into the universality class of the (supercritical) Hopf bi- 
furcation. However, when all nonlinearities are taken into 
account, the RE (|4]) give rise to heteroclinic orbits instead of 
limit cycles. The latter rapidly approach the boundaries of 
the phase space, and thus are in general well separated from 
the limit cycles predicted by (18). When comparing results 



inferred from the CGLE and stochastic lattice simulations 
in the results section, we have shown how this causes some 
quantitative mismatch, stemming from the differences be- 
tween the solutions of ([4]) and (18). However, we have also 
seen that most features of the system are actually aptly cap- 
tured by the normal form (13). Elsewhere, it will be shown 



that mvitations between subpopulations lead to limit cycles 
resulting from a Hopf bifurcation. 

3.4 Spatial structure and the continuum 
limit 

The system under consideration possesses spatial degrees of 
freedom, which are neither taken into account in the RE (|4]) 
nor in the normal form ( 13 1. Here, within a proper contin- 



uum limit, we show how the spatial arrangement of mobile 
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individuals may be included into the analytical description 
by supplementing the RE ( 4| with spatially dependent den- 
sities and diffusive terms (23|. When one additionally ac- 
counts for stochastic effects (internal noise, see below), the 
system is aptly described by a set of stochastic partial dif- 



ferential equations (SPDE) (|30l. 

The reactions (jlj and the exchange processes take place 
on a c?-dimensional hypercubic lattice (with periodic bound- 
ary conditions) of linear size L and comprising N — L"^ 
sites. The coordination number z = 2d gives the number 
of nearest neighbors of each lattice site. We set the length 
of the lattice to unity, such that the distance between two 
nearest neighbors is Sr — N~^^'^. The density of subpopu- 
lations A, B and C at time t and site r = (ri, r^;) is de- 
noted a(r,i), b{r,t), and c{r,t), respectively. According to 
the "bimolecular" reactions ([ij , the equations of motion of 
these quantities only involve neighbors located at r ± (5r • e^, 
where {ei}i=i.. is the basis of the lattice. Here, we first 
ignore all forms of correlation and fluctuations of the local 
density. While noise will be taken into account below, due 
to diffusion, correlations between neighboring sites vanish 
in the continuum limit of large systems (see below). We 
obtain the following equation for the time evolution of its 
mean (average) value a(r,t): 



dta{r,t) 



1 

- J2 {2e[a(r±<Sr-e„i)-a(r,i)] 



-t- iia{r ± 6r ■ e^, t) [l — a{r, t) — b{r, t) — c{r, t)] 
-ac{r±6r-e„t)a{r,t)Y (19) 

For an analytical description of the lattice model, fruitful 
insights are gained by considering a continuum limit with 
N ^ oo. As the lattice size is kept fixed to 1, in this limit 
the distance Sr between two neighboring sites approaches 
zero, i.e. Sr = N~^/'^ — + 0. This allows to treat r as a 
continuous variable and the following expansion is justified: 

a(r ± Sr ■ Ci, t) = a{r, t) ± Srdia{r, t) 

+ ]^Sr^dla{r,t) + o{Sr^). (20) 
With this expression, the first term on the right-hand-side 



(RHS) of Eq. ( 19 1 becomes (up to second order) 



(2e/z)^ [a{r±Sr-e,,t)~a{r,t)\ = {e/d)Sr^d^a{r,t). 
± 

If we rescale the exchange rate e with the system size N 
according to 

e^DdN'^/'^, (21) 

with a fixed (diffusion) constant D, we note that eSr^ — Dd. 
For the other terms on RHS of Eq. ( [T9| ), only the zeroth- 
order contributions in the expansion of a{r ± Sr ■ ei,t) do 
not vanish when — > oo (i.e. Sr 0). In this continuum 
limit, Eq. il9h thus turns into 



wit the local density p{r, t) = a(r, t) + b{r, t) + c(r, t). The 
equations of motion for b{r,t) and c{r,t) are obtained sim- 
ilarly. We therefore obtain the following set of partial dif- 
ferential equations (PDE): 

dta{r,t) = DAa{r,t) + ^a(r,t)[l - p{r,t)] - aa{r,t)c{r,t) , 

dtb{r, t) = DAb{r, t) + fib{r, t) [I - p{r, t)] - ab{r, t)a{r, t) , 

dtc{r,t) = DAc{r,t) + pc{r,t)[l - p{r,t)] - ab{r,t)c{r,t) . 

(23) 

The difference to the rate equations Q lies in the spatial 
dependence through diffusive terms, proportional to the dif- 
fusion constant D. 



3.5 Noise 

The discrete character of the individuals involved in the 
May-Leonard reactions ([!]) and the exchange processes are 
responsible for intrinsic stochasticity arising in the system. 
For the treatment of this internal noise, we note a time- 
scale separation between the reactions ([T|) and the exchange 
events. Namely, in the continuum limit N ^ oo, accord- 
ing to ([2]) one has e — s- 00. This means that exchanges 
occur on a much faster time-scale than the reactions {[iJ. 
Consequently, a large number of exchange events occurs be- 
tween two reactions and can be treated deterministically. 
As shown below, the fluctuations associated with the ex- 
change processes vanish as 1/A^ (for A^ ^ 00), while those 
stemming from ^ scale as 1/\/n. The latter are there- 
fore the dominating source of noise, while the former can 
be neglected in the continuum limit. To establish these re- 
sults, we consider large system sizes A^ where a stochastic 
description in terms of Fokker-Planck equations is generally 
appropriate (van Kampen 1981 Gardiner 1983). The lat- 



ter can be obtained from Kramers-Moyal expansion (i.e. a 
system-size expansion) of the underlying master equation 
(Tauber 20081, see Appendix D. In Fokker-Planck equa- 



dta{r, t) ^DAa{r, t) + pa{r, t) [l - p{r, t)] - aa{r, t)c{r, t) . 

(22) 



tions, fluctuations are encoded in a noise matrix denoted 
B. Equivalently, a set of Ito stochastic (partial) differen- 
tial equations (often referred to as Langevin equations) can 
be systematically derived. For these SPDE, the noise, of- 
ten white, is encoded in the "square root" of the matrix 
B. Namely, in this framework, the strength of fluctuations 
and the correlations are given by a matrix C, defined as 
CC^ = B. Below, we derive the relevant contributions to 
the noise matrices B and C, which lead to the appropri- 
ate stochastic partial differential equations (SPDE) of the 
system. 

Following dGardiner 1983) (Chapter 8.5), we first show 
that fluctuations stemming from pair exchanges scale as 
1/A^. Consider two nearest neighboring lattice sites r and 
r' . The rate for an individual A to hop from r to r' is 
given by ez~^a{r)\l — a{r')] (for simplicity, we drop the 
time-dependence). Together with the reverse process, i.e. 
hopping from site r' to r, this yields the non-diagonal part 
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ofB{r,r') (see e.g. ( |Tauber[ [20081 )) 



B{r, r' ^ r) 



Nz 



{a(r)[l-a(r')]+a(r')[l-a(r)]}. 

(24) 

Similarly, the diagonal entries of B read 
^(^'^) = ]|^ E {aW[l-a(r")] + a(r")[l-a(r)]}, 

n.n.r" 

(25) 

where the sum runs over all nearest neighbors (n.n.) r" of 
the site r. It follows from these expressions that 

n.n.r" 

X {a(r)[l - a(r")] + a{r")[\ - a{r)]] . (26) 

In the continuum limit, with 5r — s- 0, we use the fact that 
5r,r' Sr'^6{r — r') and obtain 

d 

B{r,r')^ ^5r'^ E [S{r - r') - 5{r ± 5re, - r')] 



{a(r)[l — a(r ± 6rei)] + a{r ± Srei)[l — 0(5")]} 



(27) 



As in Eq. (20 1, we expand 6{r ± SrSi — r') and a(r ± 6rei) 



to second order and observe that only quadratic terms in Sr 
do not cancel. With e = DdN'^/'^ and 5r = N^^^'^, we thus 
find: 



B{r,r') 



D 



drdr' [S{r-r')a{r){l-a{r)]. 



(28) 



The noise matrix B of the Fokker-Planck equation associ- 
ated with the exchange processes therefore scales as N~'^. 
In the corresponding SPDE, the contribution to noise of the 
exchange processes scales like N^^. 

We now consider the fluctuations stemming from May- 
Leonard reactions ([ij. A detailed discussion of the treat- 
ment of fluctuations arising from discrete reactions in lattice 
systems can be found in Chapter 8 of Ref . ( Gardiner 1983 1 . 



Following the derivation therein, one recovers noise terms of 
the same form as in the corresponding non-spatial model, 
although the densities now have spatial dependence. Thus, 
they may be found via a straightforward Kramers-Moyal 
expansion of the master equation describing the well-mixed 
system, see e.g. (Traulsen et al. 2005 2006 Reichenbachl 



et al. 2006 1. Here, we report the results and relegate de- 



The SPDE for the densities a{r,t),b{r,t),c{r,t) are thus 
given by the partial differential equations ( |23[ ) supplemented 
by the corresponding noise terms, which leads to: 

dta{r, t) = DAa{r, t) + ^la{r, t) [I - p(r, t)] 

- (Ta{r,t)c{r,t) +Ca£.a , 
dtb{r, t) = DAb{r, t) + /x6(r , t) [1 - p(r , t)] 

- crb{r, t)a{r, t) + Cb£,b , 
dtcir, t) = DAcir, t) + ticir, t) [1 - p(r , t)] 

-<jb{r,t)c{r,t)+Cc^c, (30) 

where A denotes the Laplacian operator, and the Gaus- 
sian white noise terms ^i{r,t) have a spatio-temporal de- 
pendence, with the correlations 

Ur,t)C,ir\t'))^S,,Sir^r')S{t~t'). (31) 



3.6 



Complex 
(CGLE) 



Ginzburg-Landau Equation 



The reaction terms appearing in the SPDE ( 30 1 coincide 
with those of the rate equations ( [4|. Above, we have re- 
cast the latter in the normal form (13 1. Applying the same 



transformations to the SPDE ( 30 ) yields reaction terms as 
in (13 1. However, owed to the nonlinearity of the trans- 



formation, additional nonlinear diffusive terms appear in 
the spatially-extended version of (13 1. In the following, the 
latter will be ignored. Furthermore, when discussing the 
spatio-temporal properties of the system, we have encoun- 
tered rotating spiral waves, whose velocity, wavelength and 
frequency have turned out to be unaffected by noise (in the 
continuum limit). An important consequence of this flnd- 
ing is that it is not necessary to take noise into account to 
compute such quantities. This greatly simplifles the prob- 
lem and, omitting any noise contributions, we focus on two 
coupled partial differential equations which are conveniently 
rewritten as a complex PDE in terms of the complex vari- 
able z — ZA + izB (see Appendix C): 

dtz{r, t) = DAz{r, t) + (ci - iuj)z{r, t) 

-C2{l-iC3)\z{r,t)\^z{r,t). (32) 

Here, we recognize the celebrated complex Ginzburg- 
Landau equation (CGLE), whose properties have been ex- 
tensively studied (Cross and Hohenberg 1993 Aranson and 



tails of the derivation to Appendix D. As the reactions ([!]) [Kramer 20021. In particular, it is known that in two dimen- 



decouple birth from death processes, the noise matrices B 
and C are diagonal. In particular, the diagonal parts of C 
read 



Ca = 



Cb = 



1 



^a{r,t) p{r,t)) + (Tc{r,t)\ 

^b{r,t) [/i(l-p(r,t)) + (7a(r,t)] 
^c{r,t)[fi{l~ p{r,t)) + ab{r,t)\ 



sions the latter gives rise to a broad range of coherent struc- 
tures, including spiral waves whose velocity, wavelength and 
frequency can be computed analytically. 

For the system under consideration, we can check that the 



CGLE ( 32 1 predicts the emergence of spiral waves which are 
stable against frequency modulation: no Benjamin-Feir or 
Eckhaus instabilities occur. As a consequence, one expects 
no intermittencies or spatio-temporal chaos, but only rotat- 
(29) ™§ spirals (Aranson and Kramer 20021. In our discussion, 



we have verifled the validity of this prediction for various 
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sets of the parameters, both for stochastic lattice simula- 



tions and for the solutions of the SPDE ( 30 1 . In particular 



the expression ( 32 1 with the parameters ci , C2 and C3 given 



by ( 15 )-( 17 1 and the properties of the CGLE allow us to ver- 



ify that the system is always characterized by spiral waves 
in the coexistence phase. Therefore, while specific proper- 
ties (size, frequency and velocity) of the emerging spirals 
depend on the values a and /i, the general (qualitative) be- 
havior of the system is already captured by the choice of 
parameters a — fi — 1. 

In the following, we describe how characteristic proper- 
ties of the spiral waves, such as the spreading velocity, the 
selected frequency and wavelength can be inferred from the 



above CGLE (32) 



3.7 The linear spreading velocity 

We have found [see, e.g., snapshots (Fig.|5])] that in the long- 
time regime the system exhibits traveling waves. Namely, in 
the steady state, regions with nearly only A individuals are 
invaded by a front of C individuals, which is taken over by 
B in turn, and so on. The theory of front propagation into 
unstable states (see, e.g., (van Saarloos, 20031 and refer- 



ences therein) is useful to study analytically the related dy- 
namics. Indeed, to determine the spreading velocity of the 



propagating fronts one linearizes the CGLE ( 32 1 around the 
coexistence state z — [i.e. the reactive fixed point of ([4|], 
which yields 

dtz{r, t) = DAz{r, t) + (ci - w)z(r, t) + 0(2^) . (33) 
It is then useful to perform a Fourier transformation: 

/•OO 

S{k,n)= drdt z{r,t)e-''' ''^'^\ (34) 



which, together with (33 1, gives the following dispersion re- 
lation: 

n{k) =uj + i{ci -Dk"^) , (35) 

where k = |fc|. As lmD.{k) > for k'^ < ci/D, the state z = 
is linearly unstable in this range of wavevectors k. This 
confirms the analysis of the spatially homogeneous RE (|4|, 
where we already found that the coexistence fixed point 
is unstable. As for other systems characterized by fronts 



propagating into unstable states (van Saarloos 20031, from 



Eq. ( 33 ) one can now compute the linear spreading velocity, 
i.e. the speed u* at which fronts (e.g. generated by local 
perturbations around z = 0) propagate. Following a classic 
treatment, whose general derivation can, e.g., be found in 
( van Saarloos 2003 1 , the spreading velocity is obtained by 
determining a wavevector k.^, according to 



dn{k) 



dk 



ImSl(fc,) _ 
Imfc* 



(36) 



The first equality singles out A:* and the second defines the 
linear spreading velocity v* . Here, we find: 



Re/c* = 



Im/c* 



The comparison of this analytical prediction with numeri- 
cal results has revealed a good agreement, as illustrated in 
Fig. [9] and discussed in the corresponding previous section. 

3.8 Wavelength and frequency 

To determine analytically the wavelength A and the fre- 
quency f2 of the spiral waves, the (cubic) nonlinear terms 
of the CGLE (32) have to be taken into account. From the 



understanding gained in the previous sections, we make a 
traveling-wave ansatz z{r,t) = ^'e^*^*^*'' '' leading to the 
following dispersion relation (with q = \q\) 



n{q) 



• + i{ci - Dq^) - C2ii + C3)Z^ 



(38) 



Separating real and imaginary parts, we can solve for Z, 
resulting in = (ci — Dq^)/c2- As already found above, 
the range of wavevectors that yield traveling wave solutions 
is therefore given hy q < \fc\fD. The dispersion relation 
can, upon eliminating Z , be rewritten as 



^{q) ^ + c-i{Dq^ ^ Cx). 



(39) 



As manifests on the RHS of ( 39 1 , comprises two contri- 
butions. On the one hand there is cj, acting as a "back- 
ground frequency" , which stems from the nonlinear na- 
ture of the dynamics and is already accounted by ([4]) when 
the system is spatially homogeneous. On the other hand. 



the second contribution on the RHS of ( 39 1 is due to the 



spatially-extended character of the model and to the fact 
that traveling fronts propagate with velocity w*, therefore 
generating oscillations with a frequency of v* q. Both con- 
tributions superpose and, to sustain a velocity v* , the dy- 
namics selects a waven umber g'^"' according to the relation 

Solving this 
< \J cxjD yields 



„scl 



van Saarloos 



ri(g>5<=') = 0^-1- v^q 

equation for (f^^ under the restriction q 



20031 

scl 



9' 



scl 



\/ci 



1 



l+d 



(40) 



Analytical expressions of the frequency ri((j'^°') and of the 
wavelength of the spirals, A = 27r/(7^°', can be obtained im- 



mediately from (39 1 and (40). In fact, the frequency reads 



scl\ 



2ci 

;+ — 
C3 

and the wavelength is given by 



27rc3Vi? 



51(1 - yrr^) 



(41) 



(42) 



The expressions (40 1-(|42| have been derived by consid- 
ering a traveling wave ansatz as described above. The 
latter hold in arbitrary dimensions. However, while trav- 
eling waves appear in one dimensions, in higher dimen- 
sions, the generic emerging structures are somewhat dif- 
ferent. E.g. rotating spirals arise in two dimensions, as 
described in this article, while scroll waves are robust solu- 



(37) tions of the CGLE (32 1 in three spatial dimensions ( Aranson 
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and Kramer 20021. However, the characteristic properties 



of these patterns, such as wavelength and frequency, still 
agree with those of traveling waves. Indeed, concerning the 
dynamical system investigated in this article, we have shown 
how the self-forming spirals are well characterized by the ex- 



pressions (41 1 and (42). The same system studied in three 



dimensions is therefore expected to exhibit an entanglement 
of scroll waves, whose wavelengths and frequencies are again 



given by Eqs. (41 1 and (42) 



4 Discussion 

Individuals' mobility as well as intrinsic noise have crucial 
influence on the self- formation of spatial patterns. We have 
quantified their influence by investigating a stochastic spa- 
tial model of mobile individuals experiencing cyclic dom- 
inance via interactions of 'rock-paper-scissors' type. We 
have demonstrated that individuals' mobility has drastic 
effects on the emergence of spatio-temporal patterns. Low 
exchange rate between neighboring individuals leads to the 
formation of small and irregular patterns. In this case co- 
existence of all subpopulations is preserved and the ensuing 
patterns are mainly determined by stochastic effects. On 
the other hand, in two dimensions, larger exchange rates 
(yet of same order as the reaction rates) yield the formation 
of (relatively) regular spiral waves whose rotational nature 
is reminiscent of the cyclic and out-of-equilibrium ensuing 
kinetics. In fact, the three subpopulations endlessly, and 
in turn, hunt each other. The location and density of the 
spirals' vortices is either determined by initial spatial inho- 
mogeneities, if these take pronounced shape, or by stochas- 
ticity. In the latter case, internal noise leads to an entangle- 
ment of many small spirals and a universal vortex density 
of about 0.5 per square wavelength. Increasing the diffusion 
rate (i.e. individuals' mobility), the typical size of the spi- 
ral waves rises, up to a critical value. When that threshold 
is reached, the spiral patterns outgrow the two-dimensional 
system and there is only one surviving subpopulation cov- 



ering uniformly the system ( Reichenbach et al. 2007a I 



The language of interacting particles enabled us to de- 
vise a proper treatment of the stochastic spatially-extended 
system and to reach a comprehensive understanding of the 
resulting out-of-equilibrium and nonlinear phenomena. In 
particular, we have shown how spatio-temporal properties 
of the system can be aptly described in terms of stochas- 
tic partial differential equations (SPDE) and confirmed our 
findings with lattice simulations. We have paid special at- 
tention to analyze the wavelength and frequency of the spi- 
ral waves, as well as the velocity of the propagating fronts. 
Numerical solutions of the SPDE have been shown to share 
(statistically) the same steady states as the lattice simula- 
tions, with the emerging spiral waves characterized in both 
cases the same wavelength, overall sizes and frequency. We 
have also studied the infiuence of stochasticity on the prop- 
erties of the coexistence state and its spatio-temporal struc- 
ture. Namely, we have compared the results obtained from 



the SPDE with those of the deterministic PDE (obtained 
by dropping the noise contributions in the SPDE), which 
still yield spiralling structures. This allowed us to shed 
light on the fact that, in the presence of (sufficient) mo- 
bility, the wavelength and frequency of the spirals are not 
affected by internal noise. However, there are major differ- 
ences between the stochastic and deterministic descriptions 
of the system. One of the most important is the infiuence 
of the initial conditions. On the one hand, if initial spatial 
inhomogeneities are larger than the noise level, or if noise 
is absent as in the deterministic descriptions, these initial 
spatial structures determine the position of the spirals' vor- 
tices. In this situation, the system "memorizes" its initial 
state, and the latter crucially infiuences the overall size of 
the emerging spiral waves. On the other hand, for rather 
homogeneous initial densities (at values of the unstable reac- 
tive fixed point), the patterns emerging from the stochastic 
descriptions (lattice simulations and SPDE) are caused by 
noise and characterized by a universal density of 0.5 spiral 
vortices per square wavelength. While we have provided 
qualitative explanations of these findings, a more profound 
understanding is still desirable and could motivate further 
investigations. 

We have also shown that analytical expressions for the 
spirals' wavelength and frequency can be determined by 
means of a complex Ginzburg-Landau equation (CGLE) ob- 
tained by recasting the PDE of the system, restricted onto 
an invariant manifold, in a normal form. There is good 
agreement between analytical predictions stemming from 
the system's CGLE and the numerical results (obtained 
from stochastic lattice simulations as well as the numerical 
solution of the SPDE). This can be traced back to the fact 
that May-Leonard rate equations are characterized by het- 
eroclinic orbits very much reminiscent of limit cycles result- 
ing from a Hopf bifurcation. The fact that the dynamics can 
be recast in the form of a CGLE, known to give rise to the 
emergence of coherent structures, reveals the generality of 
the phenomena discussed in this work and greatly facilitates 
their quantitative analysis. In particular, the emergence of 
an entanglement of spiral waves in the coexistence state, 
the dependence of spirals' size on the diffusion rate, and 
the existence of a critical value of the diffusion above which 
coexistence is lost are robust phenomena. This means that 
they do not depend on the details of the underlying spatial 
structure: While, for specificity, we have (mostly) consid- 
ered square lattices, other two-dimensional topologies (e.g. 
hexagonal or other lattices) will lead to the same phenom- 
ena, too. Also the details of the cyclic competition have no 
qualitative infiuence, as long as the underlying rate equa- 
tions exhibit an unstable coexistence fixed point and can 
be recast in the universality class of the Hopf bifurcations. 
We still note that instead of defining the model in terms 
of chemical reactions, as done here ([TJ, we can equivalently 
choose a formulation in terms of payoff matrices ( Maynard 
Smith! |1982]|Hofbauer and Sigmund) |1998[ ). 



We have investigated the system's behavior in two spa- 
tial dimensions. However, our approach, using a continuum 



16 



limit to derive the SPDE (30 1 as well as the CGLE (32 1, is 



equally valid in other dimensions and expected to describe 
the formation of spatial patterns, as long as the mobility is 



below a certain threshold value ( Reichenbach et al. 2007a I 



As examples, in one dimension, the CGLE yields traveling 
waves, while "scroll waves", i.e. vortex filaments, result in 



three dimensions ( Aranson and Kramer 2002 ) 



In this article, we have mainly focused on the situation 
where the exchange rate between individuals is sufficiently 
high, which leads to the emergence of regular spirals in two 
dimensions. However, when the exchange rate is low (or 
vanishes), we have seen that stochasticity strongly affects 
the structure of the ensuing spatial patterns. In this case, 
the (continuum) description in terms of SPDE breaks down. 
In this situation, the quantitative analysis of the spatio- 
temporal properties of interacting particle systems requires 
the development of other analytical methods, e.g. relying 
on field theoretic techniques ( Mobilia et al. 2007 ) . Fruitful 



insights into this regime have already been gained by pair 



approximations or larger-cluster approximations (Tainaka 



1994 



20071 



|Sato et al.|[T997l|Szab6 et al.jpOMl |Szab6 andTath 



The authors of these studies investigated a set of 
coupled nonlinear differential equations for the time evolu- 
tion of the probability to find a cluster of certain size in 
a particular state. While such an approximation improves 
when large clusters are considered, unfortunately the effort 
for solving their coupled equations of motion also drastically 
increases with the size of the clusters. In addition, the use 
of those cluster mean-field approaches becomes problematic 
in the proximity of phase transitions (near an extinction 
threshold) where the correlation length diverges. Investiga- 
tions along these lines represent a major future challenge in 
the multidisciplinary field of complexity science. 



Acknowledgments 

Financial support of the German Excellence Initiative via 
the program "Nanosystems Initiative Munich" and the Ger- 
man Research Foundation via the SFB TR12 "Symmetries 
and Universalities in Mesoscopic Systems" is gratefully ac- 
knowledged. M. M. is grateful to the Alexander von Hum- 
boldt Foundation and to the Swiss National Science Foun- 
dation for support through the grants IV-SCZ/ 11 19205 and 
PA002- 119487, respectively. 



from y- to z-variables yielding the normal form is detailed 
in a third appendix. 

In Subsection |3.2| we have introduced proper coordinates 
(ua, yB,yc) by y = Sx with the matrix S given by 



S 




(43) 



Hereby, the vector x = {a — a* , b — b* , c ~ c*)'^ encodes the 
deviation of the densities from the reactive fixed point. Of 
course, the time evolution of x is just given by the time evo- 
lution of the densities, Eqs. (|4|: dtX — {dta,dtb,dtc)'^ ■ For 
the temporal evolution of the y-variables, we have to apply 
the transformation given by S: dty ~ SdtX. Expressing 
them in terms of y-variables, we eventually obtain 



dtVA 



dtVB 



dtyc 



IJjCF r /— -1 r 2 2 1 

^(^^[2/A + V3yB] + ^a[y^-yB] 
- '^VaVb - ^yc [(6Ai + cr)yA - VScrys] , 

/iCr 

2(3^ + 0-) 
V3 



[ys - VSy^] 

1 



[y\~yl] 



2 o-J/Ays - ^yc [VSayA + (6// -t- cr)yB] , 
f^yc - (3m + + 7 [^A + y|] ■ 



(44) 



Using the invariant manifold, Eq. ( 12 1, we eliminate yc from 



the above and are left with equations for yA,yB alone. Ac- 



cording to Eq. ( 12 ) yc has been determined to second order 



in yA,yB, and as yc contributes to the time-evolution of 
yA, yu through quadratic terms, we obtain dtyA, dtys up to 
third order: 



dty A = 



ficr 



2(3/i - 
cr(3/i + (T 



■[yA 



[yl 



yl] 



VAyB 



8/i(3/i 
dtys = 



-2cr) 

/id 



{yl 



2(3// + a) 
a(3/i -I- cr) 



[(6/i -I- a)yA - Vicrys] + o{y^) 
ys - VSy^l 



yl) 



[vl 



vl] 



V3 



8/i(3/i + 2a 



[V3cryA + iS>ii + a)yB\ 



f^yAVB 

o{y'). 

(45) 



Appendix A. Equations for the 
variables 



These equations describe the system's temporal evolution 
y~ on the invariant manifold. 



In this Appendix, as well as the two following, some fur Appendix B. The invariant manifold 



ther details on the derivation of the normal form ( 13 ) of the 
May-Leonard RE ([4| are given. Namely, we derive the time 
evolution of the y-variables introduced in Subsection |3.2| 
and use the invariant manifold to eliminate the stable de- 
gree of freedom. The derivation of the invariant manifold 
is presented in the next section, while the transformation 



We provide further details concerning the derivation (to sec- 
ond order) of the invariant manifold of the RE (B| given by 



Eq. (121 



To determine the invariant manifold parameterized by 
yc — G{yA,yB) up to second order in yA,yB, we make 
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the ansatz G{yA,yB) — K{y\ + y^) + o{y'^), and calculate 
K. The condition (111 turns into 



2KyAdtyA + ^KysBtyB = dtyc 



VC=G 



(46) 



Using Eqs. (44 1, we obtain (to second order in yAiys) 



K- 



fMJ 



3/1 + (J 
which is satisfied for 



{yA + yB)^{-^^K+^){yl + y%) , (47) 



K 



a 3/1 + (7 
4/t 3/i + 2(7 



(48) 



Appendix C. Normal Form: the non- 
linear transformation y ^ z 



The normal form (13) of the May-Leonard RE ^ follows 
from the time-evolution equations in the y-variables, given 



by ( 45 1 , through an additional nonlinear variable transfor- 
mation. Here, we present the latter. 



The equations of motion ( 45 1 comprise quadratic and cu 



bic terms. To recast Eqs. (45 1 in their normal form, we seek 



a transformation allowing to eliminate the quadratic terms. 
We make the ansatz of a quadratic transformation y — > z 
and determine the coefficients by cancelling the quadratic 
contributions to the RE in the z variables, this leads to 



ZA 



ZB 



VA 
VB 



3/t + (T 

28/i ' 
3/i -I- a 
28/i 



2V5yAyB - 5y%] 



(49) 



To second order, this nonlinear transformation can be in- 
verted: 



VA 



VB 



ZA — [v 3z^ -I- lOzAZB 



28/1 
(3/i + (t)2 



^z%- 



ZB 



3/i -t- cr 
28/x 
(3/x + (t)2 

14/i2 



[z\ + ZAzl] + o{z^) , 
[5z^ — 2^/?,zazb — 52^] 



[zaZb + 4] + o{z^) . 



(50) 



With these expressions, one can check that equations of 



motion (45 1 are recast in the normal form (13 1 



Appendix D. Kramers-Moyal expan- 
sion of the master equation 

For a large number N of interacting individuals, the master 
equation describing the stochastic system may be expanded 
in the system size N ^ often referred to as Kramers-Moyal ex- 
pansion (van Kampen 1981 Gardiner 1983 Tauber 20081. 



As a result, one obtains a Fokker-Planck equation which is 



equivalent to a set of Ito stochastic differential equations 
(with white noise). In Subsection |3.5[ we have shown that 
for the present stochastic spatial system, only noise terms 
stemming from the reactions ([l]) contribute. The latter may 
be derived considering the stochastic non-spatial system. 
Here, we follow this approach. Starting from the master 
equation for the stochastic well-mixed system, we detail the 
system size expansion, and derive the Fokker-Planck equa- 
tion as well as the corresponding Ito stochastic differential 
equations. 

Denote s — (a, b, c) the frequencies of the three subpop- 
ulations A, B, and C. The Master equation for the time- 
evolution of the probability P{s, t) of finding the system in 
state s at time t reads 

dtPis, t) = ^ {Pis + 5s, t)Wis + Ss^ s) 

5s 

- Pis,t)Wis^ s + Ss)} . (51) 

Hereby, yV(s s + Ss) denotes the transition probability 
from state s to the state s + Ss within one time step; sum- 
mation extends over all possible changes Ss. The relevant 
changes Ss in the densities result from the basic reactions 
([T|); as an example, concerning the change in the density of 
the subpopulation A, it reads Ssa = ^/N in the reaction 
A(Z) AA, Ssa = -l/N in the reaction CA C0, and 
zero in the remaining ones. Concerning the rates for these 
reactions, we choose the unit of time such that, on average, 
every individual reacts once per time step. The transition 
rates resulting from the reactions ([T|) then read W — Ncrac 
for the reaction CA — > C0 and W ~ N^a{l ~ a ~ b ~ c) 
for AQ) — ^ A A. Transition probabilities associated with all 
other reactions ([ij follow analogously. 



The Kramers-Moyal expansion (Tauber 2008) of the 



Master equation is an expansion in the increment Ss, which 
is proportional to N^^. Therefore, it may be understood 
as an expansion in the inverse system size N~^. To sec- 
ond order in Ss, it yields the (generic) Fokker-Planck equa- 



tion (Tauber 2008): 



dtP{s,t) = -d,[a,{s)P{s,t)] + -d,dj[B,j{s)P{s,t)] . (52) 

Hereby, the summation convention implies sums carried over 
the indices i,j G {A,B,C}. According to the Kramers- 



Moyal expansion, the quantities and Bij read (Tauber 
|2008| ) 



ai{s) — ^ (5sjyV(s s + Ss) , 

Ss 

B^j{s) ^'^SsiSsjW{s ^ s + Ss). 



(53) 



Ss 



Note that B is symmetric. As an example, we now present 
the calculation of a^(s). The relevant changes Ssa = Sa 
result from the reactions A(d AA and CA C0. The 
corresponding rates as well as the changes in the density of 
subpopulation A have been given above; together, we obtain 
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aA^s) ~ fia{l — a — b — c) — aac. The other quantities are 
computed analogously; eventually, one finds 



and 



Q^a('S) — /ia(l — a — 6 — c) — aac , 
aB{s) — — a — 6 — c) — aab , 
ac{s) = /ic(l — a — b — c) — abc , 



'Baa!*) = ^ [^1(1 — a — 6 — c) + crac] , 
Sbb(s) = A^"^ - a - & - c) + cra&] , 
Beds) = N-^ [/xc(l -a~b-c) + abc] . 



(54) 



(55) 



The well-known correspondence between Fokker-Planck 
equations and Ito calculus (Gardiner 19831 implies 



that ( 52 ) is equivalent to the following set of Ito stochastic 



differential equations: 



dta = aA + Caa^a , 
dtb = as +Cbb^b , 
dtc ^ ac + Cede ■ 



(56) 



Hereby, the denotes (uncorrelated) Gaussian white noise 
terms. T he matrix C is defined from B via the relation 
CC'^ = B (Gardiner 1983 1. As B is diagonal, we may choose 



C diagonal as well, with the square roots of the correspond- 
ing diagonal entries of B on the diagonal. We obtain the 



expressions (29 1 
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